home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / interpol.pro < prev    next >
Text File  |  1997-07-08  |  3KB  |  108 lines
  1. ; $Id: interpol.pro,v 1.4 1997/01/23 00:11:23 dave Exp $
  2. ;
  3. ; Copyright (c) 1982-1997, Research Systems, Inc.  All rights reserved.
  4. ;    Unauthorized reproduction prohibited.
  5.  
  6. FUNCTION INTERPOL, V, X, U
  7. ;+
  8. ; NAME:
  9. ;    INTERPOL
  10. ;
  11. ; PURPOSE:
  12. ;    Linearly interpolate vectors with a regular or irregular grid.
  13. ;
  14. ; CATEGORY:
  15. ;    E1 - Interpolation
  16. ;
  17. ; CALLING SEQUENCE:
  18. ;    Result = INTERPOL(V, N)     ;For regular grids.
  19. ;
  20. ;    Result = INTERPOL(V, X, U)    ;For irregular grids.
  21. ;
  22. ; INPUTS:
  23. ;    V:    The input vector can be any type except string.
  24. ;
  25. ;    For regular grids:
  26. ;    N:    The number of points in the result when both input and
  27. ;        output grids are regular.  The output grid absicissa values
  28. ;        equal FLOAT(i)/N_ELEMENTS(V), for i = 0, n-1.
  29. ;
  30. ;    Irregular grids:
  31. ;    X:    The absicissae values for V.  This vecotr must have same # of
  32. ;        elements as V.  The values MUST be monotonically ascending 
  33. ;        or descending.
  34. ;
  35. ;    U:    The absicissae values for the result.  The result will have 
  36. ;        the same number of elements as U.  U does not need to be 
  37. ;        monotonic.
  38. ;    
  39. ; OPTIONAL INPUT PARAMETERS:
  40. ;    None.
  41. ;
  42. ; OUTPUTS:
  43. ;    INTERPOL returns a floating-point vector of N points determined
  44. ;    by linearly interpolating the input vector.
  45. ;
  46. ;    If the input vector is double or complex, the result is double 
  47. ;    or complex.
  48. ;
  49. ; COMMON BLOCKS:
  50. ;    None.
  51. ;
  52. ; SIDE EFFECTS:
  53. ;    None.
  54. ;
  55. ; RESTRICTIONS:
  56. ;    None.
  57. ;
  58. ; PROCEDURE:
  59. ;    Result(i) = V(x) + (x - FIX(x)) * (V(x+1) - V(x))
  60. ;
  61. ;    where     x = i*(m-1)/(N-1) for regular grids.
  62. ;        m = # of elements in V, i=0 to N-1.
  63. ;
  64. ;    For irregular grids, x = U(i).
  65. ;        m = number of points of input vector.
  66. ;
  67. ; MODIFICATION HISTORY:
  68. ;    Written, DMS, October, 1982.
  69. ;    Modified, Rob at NCAR, February, 1991.  Made larger arrays possible 
  70. ;        and correct by using long indexes into the array instead of
  71. ;        integers.
  72. ;-
  73. ;
  74. on_error,2                      ;Return to caller if an error occurs
  75. m = N_elements(v)               ;# of input pnts
  76.  
  77. if N_params(0) eq 2 then begin    ;Regular?
  78.     r = findgen(x)*(m-1)/(x-1>1) ;Grid points in V
  79.     rl = long(r)        ;Cvt to integer
  80.     s = size(v)
  81.     if s(s(0)+1) eq 1 then dif = v[1:*]-fix(v)  $ ;V[i+1]-v[i], signed for bytes
  82.     else dif = v[1:*]-v    ;Other types are already signed
  83.     return, V[rl] + (r-rl)*dif[rl] ;interpolate
  84. endif
  85. ;
  86. if n_elements(x) ne m then $ 
  87.   stop,'INTERPOL - V and X must have same # of elements'
  88. n= n_elements(u)                ;# of output points
  89. m2=m-2                          ;last subs in v and x
  90. r= fltarr(n)+V[0]               ;floating, dbl or cmplx result
  91.  
  92. if x[1] - x[0] ge 0 then s1 = 1 else s1=-1 ;Incr or Decr X
  93. ;
  94. ix = 0L                         ;current point
  95. for i=0L,n-1 do begin           ;point loop
  96.     d = s1 * (u[i]-x[ix])    ;difference
  97.     if d eq 0. then r[i]=v[ix] else begin ;at point
  98.         if d gt 0 then while (s1*(u[i]-x[ix+1]) gt 0) and $
  99.           (ix lt m2) do ix=ix+1 else $
  100.           while (s1*(u[i]-x[ix]) lt 0) and (ix gt 0) do $
  101.           ix=ix-1
  102.         r[i] = v[ix] + (u[i]-x[ix])*(v[ix+1]-v[ix])/(x[ix+1]-x[ix])
  103.     endelse
  104. endfor
  105. return,r
  106. end
  107.  
  108.